Determination of the friction coefficient of a Brownian particle 
by molecular-dynamics simulation 



F. Ould Kaddour 1,2 and D. Levesque 1 
1 Laboratoire de Physique Theorique, Universite Paris Sud, Bdtiment 210, 91405 Orsay, France 
2 Institut de Physique, Universite de Tlemcen, BP 119 Tlemcen 13000 Algerie 

(February 1, 2008) 

Abstract 

By using the Kirkwood formula, the friction coefficient of a solvated Brown- 
ian particle is determined from the integration on time of the autocorrelation 
function of the force that the solvent exerts on this particle. Extensive molec- 
ular dynamics simulations show that above a definite size of the studied sys- 
tems the value of the integral defining the friction coefficient goes to a quasi 
constant value (a plateau) when the upper bound on time increases. The 
minimal value of the system size where the integral exhibits this asymptotic 
behavior, rises with the Brownian particle size. From the plateau, a reliable 
estimate of the friction coefficient is obtained. 

PACS numbers: 83.10.Mj, 83.10.Rs 
In solutions, it is supposed that the large particles such as micelles or colloids which 
coexist with the atoms, ions or small molecules of the solvent behave as Brownian parti- 
cles. At low concentrations of the large particles, by using multi-scale analysis |I]-f§|, this 
hypothesis has been justified in the limit where the ratio between the mass of the solvated 
particles and that of the solvent molecules goes to oo. It has then been established that the 
diffusion coefficient of Brownian particles can be computed in term of the friction coefficient 
characterizing the force exerted on them by the solvent. When the Brownian particles have 



a quasi-macroscopic size, from hydrodynamic arguments the Stokes || law can be derived. 
It gives an expression of the friction coefficient £ = CRrj where R is the size of the particles, 
r) the viscosity of the solvent and C a numerical coefficient depending on the possible choices 
of boundary conditions at the interface between solvated particles and solvent. 

However in many suspensions, such as ionic solutions, the values of ratios of masses and 
sizes between the solvated particles and solvent molecules are only of the order of 10 and 
the possibility that the solvated particles can be considered as Brownian particles becomes 
questionable. Several theoretical works PJ^] and works based on numerical simulations 
|| have been devoted to this question. The main problem addressed in these works was 
that of the determination of the lower bounds of the size and mass ratios above which, to a 
good approximation, the motion of the solvated particles is Brownian. The criterion chosen 
to locate these bounds was that the diffusion coefficient of the solvated particles obeys 
to the relation between the diffusion coefficient D and friction coefficient £ strictly valid 
only for brownian particles D = k^Tj^ (&b is Boltzmann's constant, T the temperature 
of the solvent). The main concern, when the Stokes estimate of £ is used, is the choice 
of the hydrodynamic boundary condition between solvated particles and solvent which is 
well defined only when the solvated particle has a macroscopic size. This last shortcoming 
can be overcome by computing £ from its exact expression for a Brownian particle derived 
by Kirkwood |l2j and later, more rigourously, from multi-scale analysis |2]Jl3|. Obviously 
this method seems the correct way to proceed in order to check the brownian behavior of a 



solvated particle. However, as it was discussed in the literature [^|T0[|T2|,|T4|-[T^j , this method 
is not easy to use in simulations due to important finite size effects. This work is devoted 
to discuss this problem and to establish in what conditions, in a numerical simulation, the 
friction coefficient of a solvated particle of large mass and size can be credibly determined. 

The friction coefficient £ is given in terms of the integration on time t of the equilib- 
rium autocorrelation function < F(0) . F(t) > of the instantaneous microscopic force F(i) 
experienced by the Brownian particle : 
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«=i5^f <F <°>- F « >,B - (1) 

This expression of the friction coefficient has the same form as the Green-Kubo relations 
used to calculate the transport coefficients. 

For finite size systems, the computation of £ from Eq. (1) and, more generally, that of 
the transport coefficients from Green-Kubo relations are confronted with a problem that we 
illustrate for this specific case. From the momentum conservation, the force F(t) acting on 
one Brownian particle in a solvent of N molecules is given by 

F W = -^-P( t ) = -|>^, (2) 

where m is the mass of the solvent molecules with velocities v;(t) (i = 1, N) and £ N can 
be written as 

6v = Jim 6v(t) (3) 



t^oo 



= ^T* lil £ f <F (°)- F ( T ) dT 

j ft d \ rs 

= — — lim / dr- — / F(u).P(t + u) 



- _ i im < Pft)-P(O) >^v < P(0).P(0) >n (A) 
*-oo 3k B T 3k B T ' 1 ' 

The last term of Eq. (4) vanishes by symmetry on time and it is expected that, for 
t — > oo, the first term also vanishes due to the loss of correlations between P(oo) and 
P(0), according to the ergodic postulate of the equilibrium statistical mechanics, with the 
final result that £at should be zero. The well known way to overcome this paradox is that 
the thermodynamic limit on iV must be taken before that the limit t — > oo is performed. 
The argument can be summarized by guessing that at large values of t and N, £jv(t) can 
be written in the form cg(at/N) where c and a are coefficients independent of iV and 
g(at/N) is a decreasing function of t equal to at t = oo and normalized to 1 at t — 0. 
By perfoming the thermodynamic limit iV — * oo before the limit t — > oo, £ is given by 
limbec {limjv->.oo £iv(i)} = c and £ is now finite. 
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If in simulations N is large enough, it can be expected, following the remark made by 
Kirkwood in [|TjJ, that, in the range of values of t where £jv(£) reaches its asymptotic form ~ 
cg(at/N), t is such as t « a/N. In this domain of t, £tv(^0 * s gi ven by £jv(£) — c+acg'(0)t/N 
and presents a quasi plateau or a slow linear decay with t from which the value of £ in the 
thermodynamic limit can be estimated. 

As it was proposed, for instance, in HBIJlOfl it is possible to give a specific analytic form to 



g(at/N) by supposing that, following the Onsager's principle, the regression of the fluctua- 
tions of F(t) at large t is governed by the laws of the macroscopic hydrodynamics. According 
to these laws, the force exerted by the solvent on the Brownian particle is proportional to 
the momentum P(t) carried by the solvent, i.e. 

F(«) = J^PM (5) 

a relation which implies that 

< P(t) . P(0) > N = 3Nmk B Texp(-^-t) (6) 

and then 

<P(t).P(0) > N £ 
Ut) = 3^ = ^exp(-^-t). (7) 

At large N and t such as t « ^ /Nm, an expansion of the exponential yields to a linear 
expression, similar to that of £jv(i) given above, allowing to determine the friction coefficient 
as 

6v(*) = £ (i-j^i +..-)■ (8) 

In order to investigate the possibility of a computation of £ following the procedure 
describes above, we have realized a set of molecular-dynamics simulations with increasing 
values of N. 

The studied systems are made of N molecules of solvent enclosed in a periodic cubic box 
of volume V. In this box, a particle is immersed and supposed to have a size and mass M 



large compared of those of the solvent molecules. Hence the mass of this particle satisfies 
to the condition required so that the relation between D and £, given above, applies. When 
M is large, it is possible to consider that, in the time scale accessible in a simulation, the 
particle is immobile. 

The molecules and the fixed particle interact through a Lennard- Jones (LJ) potential 



modified with a cubic spline, as describe in a previous paper ||1 1|| . This potential has the 
form Vij(r) = £ijf{r/<Jij) where the indices i,j = 1 or 2 refer to the solvent and fixed 
particle, respectively. The parameters cr^ are such as ayi = (an + cr 2 2)/2, and are equal 



£12 = e 22 = en- The unit of time is chosen equal to r = J (ma n /en) and that of energy, 
length, and mass are chosen, respectively, equal to en, an, and m. The values of the solvent 
density and temperature are p* = Na n /V ~ 0.84 and T* = k B T/en = 1.0, thus specifying 
a dense liquid state near the triple point of the LJ system. The simulations were realized 
at constant energy using the standard Verlet algorithm [pH , with a time step At = 0.005 To- 
Typical simulation runs are carried out for 20000 equilibration time steps followed by 4 to 8 
millions time steps. During the runs, the time autocorrelation function of F(t) is computed 
over a sequence of blocks of 4000At to allow an evaluation of statistical errors. We have 
considered two different sizes for the Brownian particle, namely 022 = 4.0 and 7.0 and 
systems of increasing values of N : 864, 1500, 5324, 12000, 32000 and 55296. In order to 
maintain constant the value of the pressure, the volume of the simulation box was slightly 
increased, when 022 was varied from 4.0 to 7.0. 

We first discuss the case of 022 = 4.0. In Fig. 1, we show £jv(t) as a function of reduced 
time t/r . When N is increased from 864 to 32000, the behavior of £/y(t) changes drastically 
in the domain of t/ro > 5.0. For the low values of N, Civ(^) goes rapidly to zero and for 
the two larger values of N, it is almost constant. Qualitatively, this behavior of £jv(t) at 
large time corresponds to that expected when iV increases. In particular, if this behavior 
is described by Eq. (7), the results presented in Fig. 1 can be interpreted as the transition 
between the exponential decrease of £jv(t) given by Eq. (7) and the linear decrease, at large 
N, given by Eq. (8). Quantitatively, it should be possible to obtain the value of £ from the 



fit of £jv(i), f° r V r o between 12 and 20, to an exponential form when N is ~ 1000 and a 
linear form when iV is ~ 20000. 

However, the possibility that the values of £, determined from the fits made at different 
values of N, coincide within the statistical uncertainties, supposes that finite size effects, in 
particular those associated to the use of the periodic boundary conditions, do not affect the 
long time behavior of £jv(t). From the fit of the multiplicative constant of the exponential 
(cf. Eq. (7)) the estimates of f are : 219.4 (N = 864), 144.7 (N = 1500) and 122.9 
(N = 5324) and from that of the constant term in the linear form (cf. Eq. 8) they are : 
113.3 (N = 12000) and 110.7 (N = 32000). Clearly, for the two large systems these values of 
£ agree within the statistical error equal to 10-15%. But the factor of 2 between the values 
found at N = 864 and N = 32000 must be attributed to finite size effects. The side lenghts 
L of the simulation cells being for these two values of N equal to ~ 10 and ~ 35, due to 
the periodic boundary conditions the fixed particle is distant from these nearest replica by 
the same lengths. The sound velocity c s of the solvent, for the considered thermodynamic 
state, being in reduced units ~ 6, gives typical times L/c s of 1.5 and 6 beyond of which 
the finite size effects resulting from the mutual influence between the fixed particle and its 
replicas can affect the correlation function. The magnitude of these effects on the value 
of £ is difficult to assess. It has been quantitatively discussed only for the bulk values of 
the transport coefficients of diffusion or viscosity in dense fluids [|l8j, corrections of about 
~ 10% have been found for systems of N ~ 1000 molecules and they seem much larger on 
£. Other estimates of £ are obtained from the fit of the coefficient of t in the exponential 
(Eq. (7)) or the linear approximation (Eq. (8)). For the small values N, we obtained 205.5 
(N = 864) and 143.3 (N = 1500), for the large values of N, 98.4 and 91.1. These results 
seem to confirm that the asymptotic behavior of is well described by Eq. (7) taking 

into account the finite size effects. 
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FIGURES 
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FIG. 1. £jv(i) for the fixed particle of size 022 = 4.0 and increasing values of N. Insert : 
asymptotic behavior with error bars and its linear fit at N = 12000 and 32000 
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FIG. 2. £jv(t) for the fixed particle of size G22 = 7.0 and increasing values of N. Insert : 
asymptotic behavior wih error bars and its linear fit for N = 32000 and 55296 
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We consider now the case of the fixed particle with a size er 2 2 = 7.0. The £jv(0 functions 
are plotted in Fig. 2. for N equal to 5324, 12000, 32000 and 55296. The behavior of£ N (t), 
when N increases, is similar to that obtained with the particle of size 022 = 4.0, but it is only 
for iV > 32000 that presents a slow linear decrease for t/r > 10.0. For the particle 

of size a 22 = 4.0, we remark that this latter behavior is reached for systems smaller by a 
factor ~ 2 and, then, the size of the fixed particle has an important influence on the value 
of the system size where £at(£) exhibits a slow linear decrease at large time. This remark is 
confirmed by the results of the fits of £at(0 by an exponential at N = 5324 and 12000 or a 
linear function at iV = 32000 and 55296. The values of f are 217.0, 244.2, 209.3 and 207.5 
respectively. We notice that, the linear behavior of 6v(t), at large t, corresponds obvioulsy 
to the fact that < F(0) . F(t) > decreases very slowly in the same domain of t. This point 
is illustrated in Fig. 3 where this behavior is clearly seen for iV = 12000, 32000 and 55296 
within statistical errors. 
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FIG. 3. Autocorrelation function < F(0) . F(t) >n for the particle of size 022 = 7.0 at 
N = 12000, 32000 and 56296. Solid lines : estimate of its almost constant value from its av- 
erage for t/T Q > 12.0 and < 20.0 : ~ -12.6, -6.6 and -6.1. 



From the results of the present simulations it seems needed to adopt a critical point of 
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view on the previous works made in order to determine the friction coefficient of a brownian 
particle and to check when a particle of large masses and sizes can be considered as a 
brownian particle. The most important criticism is that, in these works the sizes of the 
system studied in the simulations were too small. This remark applies, for instance, to the 
simulations presented in || where the friction coefficient of a fixed hard sphere of diameter 
Ad in a solvent made of hard spheres of diameter d was computed at a density pd 3 ~ 0.471. 
Such a computation corresponds closely to that made in this work for a fixed particle with 
(J22 = 4.0. In H, the largest considered system had a size of N = 1500 which, as discussed 
above, seems too small to obtain a good estimate of £ from an exponential fit of ^{t) at 
large times and, then, to check the validity of Stokes estimate of £. The system size used 



in |10| for the computation of the friction coefficient of a fixed particle in a LJ type system 
being of the order of N = 1000, finite size effects should also affect the simulation data. 
As mentionned already, the asymptotic form of £jv(t) at large times has been discussed in 



many works in the literature, for instance in |TJJ , |15| and jTj| , in particular the question of 
the occurence of a domain of time where the friction coefficient £j\r(£) should exhibit a plateau. 
It has been proposed in || to bypass the search of such a plateau in £jv(0 by computing 
£ from the integration of < F(0) . F(£) > from t = to the value of t = t\ where, for the 
first times when t increases, < F(0) . F(£i) > becomes 0. In our simulations ti/r ~ 0.5, 
clearly from the comparison between Fig. 2 and Fig. 3 of £jv(t) and < F(0) . F(i) > N , such a 
method to estimate £ seems problematic since the value of £jv(ii) for instance at N = 55296 
does not agree with the value £ obtained from the analysis of the asymptotic behavior of 

6v(*). 

Since the present simulations show that system sizes of N ~ 20000 are needed to correctly 
estimated £, it is expected that similar system sizes are needed to compute D in order to 
avoid finite size effects. For instance in []TJ] for iV = 5324, for a brownian particle with 
(J22 = 4.0 and M = 60 in a LJ solvent at a thermodynamic state identical to that considered 
here, it was found D = 0.077. By using D = fc^T/£ this value of D agrees well with that of 
£ ~ 120 obtained in this work at N = 5324, but it differs by 10% from that, computed at 
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N = 32000, £ ~ 110. A new determination of D at iV = 32000 gives D = 0.095, which now 
agrees with our estimate of £ at the same value of N from the slow linear decay of ^(t) for 
t/r > 12. 

In conclusion from our simulations realized for increasing system sizes from N ~ 1000 to 
~ 56000, we have shown that the qualitative behavior of < F(0) . F(t) > and, consequently, 
that of is in excellent agreement with that guessed by Kirwood fl2| and in subsequent 
works. |l5|Jl6l . For large times, the representation of £jv(t) by Eq. (7), derived from a simple 
argument based on Onsager principle, seems an adequate model of this asymptotic behavior. 
From a quantitative point of view, in spite of simulation runs totalizing 4 to 8 millions of 
time steps the statistical uncertainties on the friction coefficient stay of the odrer of 15%, a 
reduction of this uncertainty by a order of magnitude seems beyond what it is possible to 
make by using present standard computers. We stress that the aim of this work was the 
investigation of the variation of £tv(£) with N. In the thermodynamic limit, the asymptotic 
behavior of £j\r(t) should be an algebraic decay at very large time. It was not considered 
here, in particular, because its amplitude is much smaller than the present uncertainties on 
£jv(£) in this domain of time |TT| , |I8"|] . 
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